Emergence and decay of turbulence in stirred atomic Bose-Einstein condensates 
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We show that a 'weak' elliptical deformation of an atomic Bose-Einstein condensate rotating at 
close to the quadrupole instability frequency leads to turbulence with a Kolmogorov energy spec- 
trum. The turbulent state is produced by energy transfer to condensate fragments that are ejected 
by the quadrupole instability. This energy transfer is driven by breaking the two-fold rotational 
symmetry of the condensate. Subsequently, vortex-sound interactions damp the turbulent state 
leading to the crystalization of a vortex lattice. 



Two-dimensional (2D) turbulence has been explored in 
diverse areas such as soap films Q , magnetohydroynam- 
ics 0, and meteorology (see Q and references therein), 
and can often display additional features not present in 
3D Q. However, the complexities of real fluids mean 
that the theoretical predictions are often at odds with 
observed spectra. In this regard there may be some ad- 
vantages in studying superfluids where the absence of 
viscosity and the quantization of vorticity can simplify 
the theoretical picture. For example, superfluid turbu- 
lence in liquid Helium 5] is found to exhibit analogous 
features to classical turbulence, in particular a Kolomor- 
gorov energy spectrum Q . Even more amenable to theo- 
retical description are atomic Bose-Einstein condensates 
(BECs) Q . In addition, atomic BECs allow the flexibility 
of studying the transition between 2D and 3D turbulence. 

Recent experiments on atomic BECs have generated 
vortex lattices by thermodynamically condensing a ro- 
tating thermal cloud |8| and mechanical rotation of the 
condensate in an anharmonic trap |?J 0, . In the lat- 
ter case, a quadrupolar collective mode of the condensate 
is excited. A dynamical instability 0, El leads to the 
nucleation of vortices, which subsequently crystallise into 
a lattice configuration. The time-scale for vortex lattice 
formation has been shown to be insensitive to tempera- 
ture [HI [H) , suggesting that the process is a purely dy- 
namical effect. The formation of the lattice has been sim- 
ulated using the time-dependent Gross-Pitaveskii equa- 
tion (GPE) in 2D with the inclusion of damping effects 
and in 3D ^tJ- However, questions remain over 
the dynamics involved, for example, how the dynamical 
instability seeds vortex nucleation, what is the damping 
mechanism leading to the crystallization of the lattice, 
and the role of dimensionality in the process. 

In this paper, we present evidence that 2D turbulence 
is a key feature of current experiments on vortex lattice 
formation in atomic BECs. In particular we show that 
the route to lattice formation can be divided into four 
distinct stages, as illustrated in Fig. 1. These stages are: 

- Fragmentation: The quadrupolar mode breaks down, 
ejecting energetic atoms to form an outer cloud. 

- Symmetry-breaking: The two- fold rotational sym- 
metry of the system is broken in a macroscopic man- 
ner, allowing the rotation to couple to additional modes, 
thereby injecting energy rapidly into the system. 

- Turbulence: A turbulent cloud containing vortices 
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FIG. 1: The condensate energy (black line) and energy of 
the outer cloud (dashed line), defined as the region with n < 
0.2no, as a function of time, with a trap rotation frequency 
Q, = 0.77ou r - The four distinct stages of the evolution (see 
text) are indicated. 



and high energy density fluctuations (sound field) is 
formed, with a Kolmogorov energy spectrum. 

- Crystallisation: The loss of energy at short length- 
scales coupled with the vortex-sound interactions [T^, fl9| 
allow the system to relax into an ordered lattice. 

Our analysis is based on the classical field method 
whereby the GPE is used to describe both the conden- 
sate mean- field ip and fluctuations [2(J. Fluctuations in 
the initial state speed-up the evolution through the four 
stages but do not change the qualitative behaviour. For 
illustrative purposes we use an initial state without ex- 
citations. The rotating trap strongly polarizes the dy- 
namics along the z-axis, such that 2D dynamics in the 
rotating plane dominate the system. Indeed we have ver- 
ified that the same results are obtained in the flattened, 
quasi-2D geometry by solving the 3D GPE. We there- 
fore proceed in the 2D limit by solving the computation- 
ally advantageous 2D GPE, as in our previous work 
We have performed simulations for a range of condensate 
sizes, and various grid and box sizes, and find the same 
qualitative results throughout. We characterise the state 
of the system in terms of the energy, which is calculated 
using the time-independent GPE energy functional, 
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and subtract the corresponding energy of the initial state. 
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Here V ext is the external potential, m is the atomic mass, 
and g = 47rh 2 Na/m, where N is the number of atoms 
and a is the s-wave scattering length. 

To apply a rotation we use a 'weakly' elliptical poten- 
tial similar to the experiments [ll| (the stirring po- 
tential used in Ref. is strongly anharmonic). The 
potential can be expressed as, 

V ex t(x, y, t) = \mw 2 r [{x 2 +y 2 )+e {X 2 - Y 2 )} . (2) 

The first term represents the static harmonic trap with 
transverse frequency uo r . In the second term, the pa- 
rameter e characterises the trap anisotropy in the co- 
ordinate system (X, Y) which is aligned with the static 
(x, y) coordinates at t = and subsequently rotates at 
angular frequency Q. We employ a trap deformation of 
s = 0.025 9]. Whereas previous theoretical investiga- 
tions have considered the rotating frame 0, 0, ^3 ? we 
perform simulations in the laboratory frame. 

Our units of length, time and energy are the healing 
length £ = h/ ^mfi, (£/c), and the chemical potential \i — 
nog. Here c = y u/m is the Bogoliubov speed of sound, 
M = n o9 is the chemical potential and no is the peak 
condensate density. Using typical parameters for a 87 Rb 
condensate, the units of distance and time correspond to 
£ ~ 0.3 am and (£/c) ~ 10 _4 s, respectively. 

Fragmentation: At t = 0, the rotation is turned on. 
The rotating trap couples energy into the condensate by 
exciting a quadrupolar shape oscillation, while the axes 
of the quadrupole rotate at the trap rotation frequency 
Q. This excitation mimics rotation yet the system re- 
mains irrotational. In the radially-symmetric system, 
the quadrupole mode is predicted to have a resonant fre- 
quency at ft = u r /y/2 [2l[, but this is shifted slightly 
by the ellipticity. Away from the resonant frequency, the 
quadrupole oscillations have reduced amplitude and are 
stable. The condensate cycles between the initial circu- 
lar state and a higher energy elongated state, but over a 
complete cycle there is no net increase in energy. Such 
quadrupole oscillations have been observed experimen- 
tally [23. 

The condensate is predicted to be dynamically unsta- 
ble at the quadrupole resonant frequency . The insta- 
bility arises because the amplitude of the mode becomes 
so large that the quadrupolar irrotational flow can no 
longer be supported. As the condensate relaxes from the 
point of maximum elongation, the fluid cannot adjust 
sufficiently quickly back towards the centre of trap. This 
results in the shedding of fluid from the ends of the con- 
densate forming low density tails and giving the conden- 
sate a spiral shape (Fig. 2(a)(ii)). The ejected material 
collapses back onto the condensate edge, forming phase 
dislocations with the main condensate. This generates 
surface waves and 'ghost' vortices An outer, low 

density (~ O.lno) cloud is formed (Fig. 2 (a) (hi)). After 
the condensate has shed material (Fig. 2(b), solid line) 
its energy no longer decreases back to the inital value 
- energy has been transferred irreversibly into the con- 




FIG. 2: Fragmentation: (a) Snapshots of the condensate den- 
sity in the range (0 - 0.05n ) at times t w (i) 400, (ii) 800, 
(iii) 1200, and (iv) 6000 (£/c). Each plot represents a re- 
gion [—30, 30]£ x [—30, 30]£ (while the numerical box is much 
larger). Dark represents high density, (b) Total condensate 
energy (solid line), energy of inner cloud (dotted line), and 
energy of outer cloud (dashed line). 



densate. We monitor the relative evolution of the outer, 
low density cloud and the inner, high density conden- 
sate by defining the outer cloud to be where the den- 
sity is less than a certain value, taken here to be 0.2no- 
The injected energy goes primarily into the outer cloud 
(Fig. 2(b), dashed line). Note that, although the energy 
of the inner and outer clouds are comparable, the outer 
cloud contains only around 10% of the total atoms, and 
therefore the average energy per atom is considerably 
higher. Subsequently, the inner cloud continues to un- 
dergo quadrupole oscillations and eject small fragments, 
as indicated by the energy curves shown in Fig. 2(b). 
Also, the outer cloud develops more structure at short 
length-scales (Fig. 2(a)(iv)). 

For the parameters employed here, we observe the frag- 
mentation of the cloud (and ultimately the formation of 
a vortex lattice) for rotation frequencies in the range 
0.72 < Q/u r < 0.78, which is in reasonable agreement 
with the experimental results of Madison et al. 9] . Out- 
side this range the quadrupole mode is stable, and the 
width of the unstable region increases with the trap el- 
lipticity £, which is consistent with the experimental ob- 
servations of Hodby et al. [ll| . 

Symmetry breaking: The condensate and potential 
have a two-fold rotational symmetry,which is clearly vis- 
ible in Fig. 2(a). In the experiments |y,[ll| this symmetry 
is absent. Eventually in the simulations, an asymmetry 
grows out of the numerical noise generated by modelling 
a rotation using a static square grid. We characterise this 
asymmetry in terms of an asymmetry parameter defined 
as, 

J [1^ fo, y) | 2 - |^ (-x, -y) | 2 ] dxdy 
j \i/>(x,y) \ 2 dxdy 

This asymmetry parameter (Fig. 3(b), solid line) grows 
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FIG. 3: Symmetry breaking: Evolution of (a) total conden- 
sate energy and (b) condensate asymmetry parameter for the 
case of no trap jitter (solid line), and trap jitters of 7 = 0.0001 
(dashed line) and 7 = 0.1 (dotted), (c)-(d) Snapshots of the 
total condensate density and phase, respectively, during the 
symmetry breaking stage at t = 9600 (£/c) in the absence of 
trap jitter. 



exponentially over time. When it reaches macroscopic 
levels (a ~ 0.1) we enter a symmetry-breaking phase at 
t ~ 8000 (£/c). From then on, additional modes can 
be excited, and energy is rapidly coupled into the sys- 
tem (Fig. 3(a), solid line). This energy predominantly 
excites the outer, low density cloud, as shown in Fig. 1 
(dashed line). The condensate density and phase during 
this stage are shown in Fig. 3(c)-(d). The original two- 
fold rotational symmetry is now clearly broken. Towards 
the end of this stage the outer cloud strongly couples to, 
and merges with, the inner condensate. Energy and an- 
gular momentum become transferred to the inner cloud, 
the quadrupole mode finally breaks down, and vortices 
become nucleated in the higher density regions. 

In experiments, symmetry-breaking will occur due to, 
for example, the thermal cloud, quantum fluctuations, 
trap imperfections, and fluctuating fields. We model the 
effect by allowing the trap centre to randomly jump, or 
jitter, within a region [—7, +7] x [—7, +7]. Fig. 3(a)-(b) 
show the results for trap jitters of 7 = 0.0001 (dashed 
line) and 0.1 (dotted line). We observe the same qual- 
itative behaviour as in the absence of the jitter, al- 
though the macroscopic symmetry-breaking occurs ear- 
lier when jitter is added. Even for an extremely small jit- 
ter (7 = 0.0001£, 3 orders of magnitude smaller than the 
grid size) the effect is significant, thereby demonstrating 
the importance of symmetry. Symmetry-breaking allows 
vortices to enter the condensate one by one, rather than 
in opposing pairs 15] which reduces the threshold energy 
for vortex nucleation. 

Turbulence: Following the injection of energy into the 
condensate during the symmetry-breaking stage, a highly 
excited and energetic condensate containing randomly- 
positioned vortices is formed, as shown in Fig. 4(a)- (b). 
We analyse this stage of the evolution by calculating the 
energy spectrum, shown in Fig. 4(c). During the tur- 
bulent phase the system closely follows a Kolmogorov 




FIG. 4: Turbulence: Snapshots of the turbulent (a) conden- 
sate density and (b) phase at t = 11000 (£/c). The vortices 
are characterised by a node in the density and an azimuthal 
phase change of 2tt. (b) Energy spectrum in k- space during: 
(i) the turbulent stage at t = 11000 (f/c) (bold line); (ii) 
fragmentation at t = 8000 (£/c) (intermediate grey line); and 
crystallisation at t = 20000 (£/c) (light grey line). The tur- 
bulent Kolmorogov behaviour Ek oc k~ 5 ^ 3 is indicated (dot- 
dashed line). 



energy spectrum E(k) oc k 5 / 3 over a range of A;- values, 
as shown in Fig 4(c) (bold line) for t = 11000(f /c). Such 
behaviour is a key signature of classical turbulence and 
also occurs in models of superfluid turbulence 6]. The 
departure from a k~ 5 ^ 3 — law occurs at an upper bound of 
k « 2(27r/£), corresponding to the characteristic length- 
scale of vortex-sound interactions , and a lower bound 
of k ~ 0.3(27r/£), corresponding to the size of the con- 
densate. An additional feature of 2D classical turbulence 
is that the spectrum is predicted to show k~ 3 behaviour 
following the k~ 5 ^ 3 region 4]. However, we observe a 
k~ 6 dependence in this range. The power spectrum in 
the fragmentation stage (Fig. 4(a), intermediate grey) 
and the succeeding crystallisation stage (Fig. 4(a), light 
grey) do not show a k~ 5 ^ 3 behaviour. 

Crystallisation: The turbulent state contains a dense 
sound field and vortices. In previous work we have shown 
how vortices can both radiate and absorb sound waves 
in a Bose-Einstein condensate pj| Uil- We propose 
that this vortex-sound interaction enables the randomly- 
positioned vortices in our rotating condensate to crys- 
tallise into a lattice. 

To demonstrate the transfer of energy from the sound 
field (density fluctuations) to the vortices, we divide the 
energy into a component due to the vortices Ey and 
a component due to the sound field E$. We approxi- 
mate Ey at a particular time by numerically generating 
a similar condensate (by propagating the GPE in imagi- 
nary time) with the same vortex distribution but without 
sound. The sound energy is then defined as E$ = E — Ey. 
During the fragmentation, symmetry breaking and tur- 
bulent stages, both sound energy (Fig. 5(a), dotted line) 
and vortex energy (Fig. 5(a), dashed line) are fed into 
the system. At the start of the crystallisation phase the 
sound energy is considerably higher than the vortex en- 
ergy. However as time progresses the sound energy de- 
creases with this energy being transferred to the vortices. 
Figs. 5(b)-(c) shows the condensate density at the begin- 
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FIG. 5: Crystallisation: (a) Total condensate energy (solid 
line), sound energy (dotted line) and vortex energy (dashed 
line). Snapshots of the (a) total condensate density (0 — no) 
and (b) low density (0 - 0.05n ) at times (i) 20000 and (ii) 
100000 (f/c). 



ning and towards the end of the crystallisation phase. 
We see that the conversion of sound energy into vor- 
tex energy is associated with the ordering of the vortices 
from a disordered distribution to a lattice configuration 
and a smoothing of the condensate profile (Fig. 5(b)). 
Furthermore, the outer cloud shrinks (Fig.5(c)). By the 
time t = 100000 (£/c) the vortex energy is substantially 
greater than the sound energy. At this point we observe 
a well-ordered vortex lattice. 

In our simulations, the finite grid size sets a maximum 
value in momentum space, with higher k modes having 
zero occupation [20j. We have imposed reduced values 
of the momentum cutoff and consistently observe lattice 



formation, with no marked effect on its timescale. This 
further supports the idea that the vortex lattice forma- 
tion is independent of thermal effects. 

In the vortex lattice experiments |?J [H| , the lattices 
are observed after up to Is of trap rotation. In our simu- 
lation, a noisy lattice has formed by t ~ 2s (Fig. 5(b,i)), 
while it takes several more seconds for the vortices to 
settle into a clean lattice. However, as shown in Fig. 3, 
the timescale of the fragmentation stage is sensitive to 
the degree of symmetry-breaking effects present in the 
system, as well as the trap rotation frequency and ellip- 
ticity. One would therefore expect that in a real system, 
with all its inherent imperfections, the timescale for this 
stage will be reduced. 

Note that if the rotation is terminated before the peak 
energy has been coupled into the system, the lattice still 
forms, albeit at a lower energy. This allows control over 
the number of vortices which ultimately form in the lat- 
tice. 

In summary, we have shown that 'stirring' atomic con- 
densates generates turbulence. We verify that a k~ 5 ^ 3 
power law is observed for two-dimensional superfluid tur- 
bulence. The turbulent state subsequently evolves into a 
vortex lattice by vortex-sound interactions. Future work 
will focus on turbulence in spinor systems and the effect 
of dimensionality on the turbulence spectrum. 
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